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Abstract 

•/^ \ Based on density-functional perturbation theory we have computed from first principles the pho- 

toelastic tensor of few crystalline phases of silica at normal conditions and high pressure (quartz, 
a-cristobalite, /3-cristobalite) and of models of amorphous silica (containig up to 162 atoms), ob- 
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crystalline silicon and MgO as prototypes of covalent and ionic systems. The agreement with avail- 
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^ , able experimental data is good. A phenomenological model suitable to describe the photoelastic 

' properties of different silica polymorphs is devised by fitting on the ab-initio data. 
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I. INTRODUCTION 



The dependence of the refractive index versus strain in sihca glass, namely the photoe- 
lasticity, is of interest from the fundamental point of view as well as for several technological 
applications in optics and microelectronics. For instance, photoelasticity in a-Si02 is known 
to cause a reduction of fiber Bragg gratings efficiency^ and to produce a loss of resolution of 
pure silica lenses used in photolytography^. Despite its widespread interest, a microscopic 
and quantitative modeling of photoelasticity in amorphous silica is still lacking. 

In this paper, we have computed the photoelastic tensor of crystalline and amorphous 
silica by first principles within density functional perturbation theory (DFPT)^ aiming at 
identifying the microscopic properties which rule the photoelastic response in silica systems. 
Models of amorphous silica containing up to 162 atoms have been generated by quenching 
from the melt in combined classical and Car-Parrinello molecular dynamics simulations. 
Photoelastic coefficients of quartz have been computed previously within DFPT by Detraux 
and Gonze'^, who found good agreement with experimental data. From the change of ef- 
fective charges of oxygen upon strain, the latter authors have suggested that an important 
contribution to photoelasticity in quartz comes from the dependence of oxygen polarizability 
on the deformation of the SiOSi angles upon strain. 

In this paper we have enlarged the scope of the previous work by Detraux and Gonze^ by 
computing the photoelastic coefficients of a-Si02 and of several crystalline phases (quartz, 
a-cristobalite, /?-cristobalite) at normal conditions and at high pressure from first princi- 
ples. Comparison with available experimental data and with calculations including a scissor 
correction to the electronic band gap, let us conclude that the simple local density approxi- 
mation to DFPT is suitable to well reproduce the experimental photoelastic tensor, despite 
a large error (10 %) in the dielectric constants. The ab-initio photoelastic coefficients of 
the crystalline phases of silica provided us with a database over which a phenomenological 
model of the dielectric response has been fitted. This model aided us in the interpretation of 
the ab-initio results on the photoelastic properties of amorphous silica. As described below, 
we have confirmed on a quantitative basis that the dependence of the anisotropic oxygen 
polarizability on the SiOSi angle is the key parameter for modeling the dielectric and pho- 
toelastic properties of silica. The paper is organized as follows. In section II we describe our 
computational framework. In section III we report on a preliminary test of our framework 
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on the photoelastic tensor of Si and MgO as prototypes of covalent and ionic materials. In 
section IV we report our results on the photoelastic coefficients and structural response to 
strains of several crystalline phases of silica: quartz at ambient conditions and high pressure, 
a-cristobalite, and /3-cristobalite. In section V, we first report the details of the protocols 
used to generate the models of a-Si02. Analysis of their structural and elastic properties are 
given. Then the calculated photoelastic tensor for a-Si02 is discussed and compared with 
experimental data. In section VI we present a phenomenological model of the dielectric 
properties of silica that we have fitted on the dielectric tensor of a-cristobalite at different 
densities. Its transferability was proven by comparison with the ab-initio database on the 
dielectric and photoelastic tensors of the crystalline and amorphous systems presented in 
the preceding sections. Finally, section VII is devoted to discussion and conclusions. 

II. COMPUTATIONAL DETAILS 

We have used density functional theory in the local density approximation (LDA)^, norm 
conserving pseudopotentials for oxygen^ and silicon^ and plane- waves expansion of the Kohn- 
Sham (KS) orbitals up to a kinetic cutoff of 70 Ry, if not specified otherwise. The structure 
of the crystalline phases has been determined by optimizing all the structural parameters. 
Monkhorst-Pack (MP)^ meshes have been used in the integration of the Brillouin Zone 
(BZ). The calculations on the crystalline phases have been performed by using the code 
PWSCF— . Models of a-Si02 have been generated by quenching from the melt in classical 
molecular dynamics simulations using the empirical potential by van Beest et ai.— which 
has been shown to properly describe the structural and dynamical properties of amorphous 
silica^^. We have slightly modified the interatomic potential of ref."° for short Si-0 and 
0-0 separations by adding a term of the form Vij = Aij/r^"^ — Bij/r^ with A = Q eV ■ A^^ 
and B = 20 eV ■ for Si-0 interactions and A = 2A eV ■ A^^ and B = 180 eV ■ A^ 
for 0-0 interactions. The modified potential turned out to be necessary to describe the 
high-temperature liquid (above 5000 K) where ions with opposite charges may approach 
each other very closely and fall in the well of the Coulomb potential in the absence of 
the additional repulsion here introduced. This new term does not change the structural 
properties of crystalline and amorphous phases. The simulations have been performed at 
constant volume. Ab-initio annealing of the models at 600 K for 0.75 ps has been then 
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performed by Car-Parrinelloi^ molecular dynamics simulation as implemented in the code 
CPMD^^. We have generated two amorphous models containing 81 and 162 atoms. Different 
quenching and annealing protocols have been used as discussed in section V. Integration of 
the BZ has been restricted to the F point only in the Car-Parrinello simulations. A fictitious 
electronic mass of 800 and a time step of 0.15 fs have been used. In the dynamical Car- 
Parrinello simulations we have used an ultrasoft pseudopotential^"^ for oxygen and a lower 
cutoff of 27 Ry. The structures generated by Car-Parrinello simulations have been then 
optimized with norm-conserving pseudopotentials to study the dielectric properties. We 
have computed the dielectric and photoelastic tensors within density functional perturbation 
theory^. The photoelastic tensor pij^i is defined by 

^^7/ = PijkiVki (1) 

where Eij is the optical dielectric tensor and rjki is the strain tensor. Only the electronic 
contribution is included in Sij, also indicated as e°°. Experimentally, it corresponds to the 
dielectric response measured for frequencies of the applied field much higher than lattice 
vibrational frequencies, but lower than the frequencies of the electronic transitions. The 
photoelastic coefficients have been calculated by finite differences from the dielectric tensor 
of systems with strains from -2 % to +2 %. The independent components of Pijki (21 in 
general) can be reduced by the symmetry of the system. In quartz for instance there are only 
eight independent coefficients which can be determined by applying three types of strain rjn, 
?733 and 7723. The components of the photoelastic tensor will be expressed hereafter in the 
compressed Voigt notation. 

The dielectric tensor for the different strained configurations has been calculated within 
DFPT— , as implemented in the code PWSCF and PHONONS^. However, the largest amor- 
phous model considered here which contains 162 atoms turned out intractable within the 
latter framework. Therefore, we have also computed the dielectric tensor within the Berry 
phase (BP) approach to the electronic polarizatio n^^i^^ as developed by Putrino et ai.— . In 
these latter calculations we have restricted the BZ integration to the F-point only. This 
latter technique converges to the correct DFPT result in the limit of large simulation cells 
or equivalent k-points mesh, although with a different pace with respect to standard DFPT. 
How large should the system be to reproduce the correct result within the BP approach has 
been investigated for /3-cristobalite in section IVc below. It turns out that for system size 
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of the order of 192 atoms and the F-point only, the error in the photoelastic coefficients of 
/9-cristobahte are of the same order of the errors introduced by the LDA. Photoelasticity 
of the largest model of a-Si02 in section V has been addressed within the BP approach as 
implemented in the code CPMHU^. 

The exchange-correlation functionals available in literature (LDA and generalized gradi- 
ent approximation (GGA)) usually underestimate the electronic band gap and overestimate 
the electronic dielectric tensor up to 10-15 This discrepancy can be corrected semi- 
empirically by applying a self-energy correction, also referred to as a scissor correction, 
which consists of a rigid shift of the conduction bands with respect to the valence bandsti^. 
This procedure has been used successfully to reproduce the photoelasticity of Si^^, GaAs^S 
and quartz^. However, as shown in section III and IV, it turns out that even within simple 
LDA, the error in the photoelastic coefficients is smaller than what expected on the basis 
of the error in the dielectric constant itself. The scissor correction has thus been neglected 
in the calculations on amorphous silica and on crystalline phases at high pressure. The 
calculations with the scissor corrections have been performed with the code ABINIT— . 

III. TEST CASES: SILICON AND MGO 

As test cases of our theoretical framework, prior to the application to silica, we have com- 
puted the photoelastic coefficients of crystalline silicon and MgO as prototypes of covalent 
and ionic systems, respectively. 

A. Silicon 

We have computed the photoelastic tensor for the conventional cubic supercell containing 
eight atoms. The dielectric tensor has only one component and there are three photoelastic 
coefficients, pn, pi2 and P44 independent by symmetry. A 8x8x8 MP mesh over the BZ 
corresponding to 20 k-points in the irreducible wedge is used. The KS states are expanded 
in plane waves up to a kinetic cutoff of 20 Ry. The photoelastic coefficients, computed at the 
theoretical equilibrium lattice constant (a = 5.377 A), are reported in table E] and compared 
with previous results obtained by Levine et al.-^. These authors computed the photoelastic 
coefficients at the experimental lattice parameter (a = 5.429 A) within simple LDA and by 
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adding a scissor operator, corresponding to a rigid shift of 0.7 eV of the conduction band 
with respect to the valence band. Although the scissor correction is crucial to reproduce 
the dielectric constant of silicon, it produces small changes in the photoelastic tensor both 
at the experimental and at the theoretical equilibrium densities. On the other hand, the 
photoelastic coefficients are strongly dependent on the lattice constant (cfr. table H}. Better 
agreement with experiments is achieved when the theoretical lattice parameter is used. This 
conclusion is also true for the calculation of the piezoelectric tensor in other material o^^i^^ . 
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this work 


this work 


ref. 




ref^ 
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13.20 


11.38 




11.4 


11.4 


Pll - Pl2 


-0.105 


-0.110 


-0.105 


-0.115 


-0.111±0.005 


Pll + 2pi2 


-0.058 


-0.051 


-0.085 


-0.067 


-0.055±0.006 


Pll 


-0.090 


-0.090 


-0.098 


-0.099 


-0.094±0.005 


Pl2 


0.016 


0.020 


0.007 


0.016 


0.017±0.001 


Pu 


-0.052 




-0.045 


-0.049 


-0.051±0.002 



TABLE L Dielectric constant and photoelastic coefficients of crystalline silicon. Results obtained 
by simple LDA and by adding a scissor correction of 0.7 eV are compared. Previous theoretical 
results by Levine et alM are also reported. 
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B. MgO 



MgO crystallizes in a NaCl-type lattice. We have computed the photoelastic coefficients 
for a conventional cubic supercell containing eight atoms. A norm conserving pseudopoten- 
tial of Car- Von Barth type^^ including non-linear-core corrections^^ for Mg and a 8x8x8 MP 
mesh for the BZ integration have been used. The equilibrium lattice parameter has been 
obtained by fitting the equation of state with a Murnaghan function^^ which yields a = 4.197 
A {aexp = 4.2017 A) and a bulk modulus B = 185 GPa {B^xp = 160 GPa). The dielectric 
constant is, as usual, overestimated (~ 5%), but the calculated photoelastic coefficients re- 
ported in table [H] are in good agreement with experimental dat a^^i^^ . The spread in the 
experimental data is quite large, as they have not been obtained by the Brillouin scattering 
method but by a hydrostatic pressure method which is subject to larger uncertainties. 





LDA {athec 


,) Exp^ 


Exp^ 


£ 


3.16 


2.92 


3.02 


Pn 


-Pi2 -0.26 


-0.24 


-0.248 


Pii 


-0.31 


-0.30 


-0.259 


Pl2 


-0.05 


-0.08 


-0.011 


Pii 


-0.075 




-0.096 



TABLE II: Dielectric constant and photoelastic coefficients of MgO. 

IV. CRYSTALLINE SILICA 

In the next sections we report the calculated photoelastic coefficients of three crystalline 
phases of silica at normal conditions and at high pressure: a-quartz, a-cristobalite, and 
/9-cristobalite. The collected set of data have then been used to fit the phenomenological 
model of photoelasticity discussed in section VI below. 

A. a-Quartz 

The photoelastic coefficients of quartz have been recently computed within DFPT by 
Detraux and Gonze^. This phase has thus been used as a test case for our theoretical 
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framework on silica systems. The unit cell of quartz is hexagonal (space group P3i21) 
and contains nine atoms of which two are independent by symmetry in the positions Si 
(m, 0, 0) and O {x,y,z). We have used a 3x3x2 MP mesh in the BZ integration although 
a single special point is enough to achieve convergence in the structural properties^. The 
optimized lattice parameters (a, c/a) and internal structural parameters are reported in 
table mil Actually, the agreement with experiments on structural parameters is fair, but 
somehow worse than usual for the SiOSi angle and the lattice parameter. These discrepancies 
have to be mainly ascribed to the choice of the LDA for the exchange and correlation 
functional which commonly causes an underestimation of the equilibrium volume. It was 
demonstrated that GGA is necessary to reproduce the correct energy differences among 
the phases of silica^°, however the equilibrium lattice parameters and the Si-O-Si angle 
are only slightly improved for a-quartz (agreement with experimental data is still within 
2%) and Si-0 distances slightly worsen (see table 1 of ref.— for comparison). The small 
differences between our results and those reported in ref.— depend on the different choice of 
the pseudopotentials. 

The calculated dielectric constants and independent photoelastic coefficients are collected 
in table ITTTl Only three independent strains [rju, 7733, and 7723) are necessary to compute 
all the independent pij. The theoretical pij are in good agreement with previous results 
obtained by Detraux and Gonze^ by adding the scissor correction and with experiments. 
Although the scissor operator is necessary to reproduce the dielectric constant, its effect on 
the photoelastic coefficients is less important. The scissor operator has been adjusted to 2.1 
eV in ref.— to better reproduce the dielectric constants of quartz. 

To identify the main structural changes upon strain, we have computed the response of the 
structural parameters (bondlengths and angles) to strain of ±2%. Results are summarized in 
table II Vl The stiffness of the Si-0 bonds causes the structure to respond to strain mainly by 
cahnges of the bond angles. The Si-O-Si angles undergo larger changes than the tetrahedral 
0-Si-O angles due to the strongly directional character of the sp^ orbitals of silicon. 

The dielectric tensor and photoelastic coefficients have been calculated at high pressure 
as well. Table IVl reports the structural and dielectric properties of a-quartz at 3 and 7 GPa, 
well below the transition pressure to the monoclinic phase (21 GPa)'^^. The variation of 
the equilibrium volume and of the internal structural parameters of a-quartz as a function 
of pressure is in good agreement with experimental. The Si-O-Si angle undergoes larger 
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0.24 


95 


97 


F66 


0.085 


0.080 


0.11 


0.10 




0.25 


0.25 


0.27 


0.29 


Pl4 


-0.031 


-0.023 


-0.03 


-0.03 


P41 


-0.040 


-0.034 


-0.028 


-0.047 


P44 


-0.056 


-0.062 


-0.061 


-0.079 



TABLE III: Structural parameters, dielectric tensor and photoelastic coefficients of a-quartz as 
computed in this work and measured experimentall} "?^!?^!?? . Previous results by Detraux and 
Gonze^ with the scissor correction are also reported for sake of comparison. Angles are in degree 
and lengths in A. 
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0.024 
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-Si (1) 
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77.7 
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160.4 




0-Si- 


-0 (1) 


77.2 


47.5 


0-Si- 


-0 (2) 


-33.9 


-0.6 


0-Si- 


-0 (2') 


-19.0 




0-Si- 


-0 (3) 


36.0 


-51.7 


0-Si- 


-0 (3') 


-12.9 




0-Si- 


-0 (4) 


-15.7 


57.4 



TABLE IV: The response of Si-0 bondlength and Si-O-Si angles of a-quartz to r/n, r/33 and r/23 
strains. Symmetry breaking due to jyn strains generates three independent Si-O-Si angles. 

changes upon compression than the 0-Si-O angles, while the Si-0 distances are nearly 
unaffected. 

We are not aware of any experimental data on the photoelastic coefficients of quartz at 
high pressure. 

The dependence on pressure of the response of structural parameters to strain is reported 
in table IVIl By inspecton of table IVII we can conclude that at high pressure the Si-0 
bonds become even stiffer. However, since the flexibility of the Si-O-Si angle is reduced, 
also the deformation of the tetrahedral unit becomes sizable. The data reported so far 
are not sufficient to establish a correlation between structural and photoelastic properties, 
however these results support the idea that bond angles (mainly Si-O-Si) play a key role in 
determining the photoelastic coefficients of quartz. 

In order to identify the structural response of silica polymorph to strain, it is useful to 
consider this system as a network of SiOSi units. The structural response to strain can then 
be expressed in terms of the length and orientation of the vector distance Si-Si between the 
two silicon atoms of the SiOSi unit which depends on the Si-O-Si bond angle. The orientation 
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Si-O-Si 


137.7 (143.7) 
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0-Si-O 
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0-Si-O 
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£11 


2.590 


2.694 


2.812 


£33 


2.620 


2.737 


2.865 
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0.170 


0.174 


0.219 


P12 


0.227 


0.214 


0.209 


Pl3 


0.245 


0.215 


0.199 


P33 


0.085 


0.067 


0.051 


P31 


0.253 


0.251 


0.254 


P44 


-0.056 


-0.055 


-0.053 



TABLE V: Structural parameters, dielectric and photoelastic tensor of a-quartz at high pressure. 
Experimental data from Refi^ are given in parenthesis. Angles are in degree and lengths in A. 

of the SiOSi units can be valuated by computing the projection of the Si-Si distances on the 
Cartesian axis. This analysis (tab. IVIIll shows that a tensile strain increases the alignment 
of the Si-Si vectors along the strain axis. 
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TABLE VI: The response of the structural parameters of a-quartz to the 7733 strain at different 
pressures. Angles are in degree and lengths in A. 
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P:. (Si-Si) 
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1.505 
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0.000 


-0.024 


-0.006 




P.(Si-Si) 


0.000 


-0.008 
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Px(Si-Si) 
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-0.416 




Py (Si-Si) 


0.000 


0.000 


-0.250 




P.(Si-Si) 


1.762 


1.595 


1.718 



TABLE VII: The derivatives of the projection of the Si-Si vectors along the direction i (Pj(Si-Si)) 
with respect to strains rjn and 7733 in a-quartz at different pressures. 

B. a-Cristobalite 

The tt-cristobalite crystal has tetragonal symmetry (space group P4i2i2). The unit cell 
contains 12 atoms of which two are independent by symmetry at the positions Si {u, u, 0) 
and O {x,y,z). The optimized structural parameters of a-cristobalite are compared with 
experiments in table IVIIIl A 3x3x2 MP mesh has been used in the BZ integration. The 
experimental value of Si-O-Si angle (144.7°) which is suggested to be the structural param- 
eter which mostly influences the photoelastic respons^ is very similar to that of quartz 
(143.9°). Conversely, a-cristobalite and quartz differ in density (2.36 g/cm^ and 2.65 g/cm^ 
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for a-cristobalite and a-quartz, respectively) and in the ring topology (6- and 8-membered 
rings in quartz and only 6-membered rings in a-cristobalite). 
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109.36 
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110.35 
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0-Si-O (4) 
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TABLE VIII: Calculated and experimental structural parameters of a-cristobalite. Angles are in 
degree and length in A. 

The dielectric tensor has two independent components (en = £22 and £33). The indepen- 
dent photoelastic coefficients are seven and can be obtained by applying of four strains: 7711, 
^33) V12 and ri23. The theoretical dielectric and photoelastic tensors are reported in table 
IIXI The results obtained by including a scissor correction equal to that used for quartz 
(2.1 eV)^ are also reported. The scissor correction largely improves the dielectric constants, 
but produces smaller changes in the photoelastic coefficients as occurs in quartz. Hereafter, 
the scissor correction will thus be neglected in the calculation of the photoelastic tensor for 
other crystalline and amorphous phases of Si02. 

We are not aware of any experimental measurement of the photoelastic coefficients in 
a-cristobalite. 

The differences in density and topology of quartz and a-cristobalite determine a different 
response of the structural parameters to strain (results are displayed in tab. |X|). In a- 
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LDA LDA+scissor 



£11 2.363 2.206 (2.211) 



£33 2.358 2.196 (2.202) 



pn 0.218 



0.225 



Pi2 0.244 



0.248 



Pi3 0.293 



0.299 



P33 0.152 



0.162 



P31 0.294 



0.298 



P44 -0.094 



-0.090 



P66 -0.068 



-0.066 



TABLE IX: Dielectric and photoelastic tensors of a-cristobalite within the simple LDA and LDA 
supplemented by the scissor correction (2.1 eV). Experimental values for £ii and £33 are reported 
in parenthesis^. 

cristobalite bond lengths are affected by strain even less than in quartz and the distortion 
of the tetrahedral bond angles is smaller. In fact, a-cristobalite has a more open structure 
than quartz and it is able to accommodate strains mainly by a rotation of the tetrahedra. 

The derivatives of the projection of the Si-Si vectors along the Cartesian axes with respect 
to strain are reported in Tab lXIl The alignment of the Si-O-Si units along the strain axis 
observed in a-quartz occurs also in a-cristobalite for strain along the main symmetry axis 
(^33)) but it is much smaller for the rjn strain. 

C. /3-Cristobalite 

The /3-cristobalite phase can be obtained by heating a-cristobalite above 270 °C at nor- 
mal pressure^. The space group Fd3m assigned experimentally to /5-cristobalite has been 
interpreted in the past as an average structure composed of small domains with six possible 
orientations of the structure with M2d symmetryi^. More recent NMR data have shown 
that the coexistence of different domains of lower symmetry is in fact dynamical and not 
static as due to domains with fixed orientation^. However, the /3-cristobalite structure with 
/42(i symmetry is locally stable at low temperature (0-300 K) in the small simulation cells 
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Si-6-Si (1) 


-31.2 


70.6 


±33.3 


Si-6-Si (2) 


186.7 




±37.2 


0-Si-O (1) 


9.98 


-35.3 


±30.5 


0-Si-O (1') 


68.2 




±18.2 


0-Si-O (2) 


-14.9 


3.81 


±41.9 


0-Si-O (3) 


-16.9 


14.5 


±20.5 


0-Si-O (3') 


-18.4 




±54.4 


0-Si-O (4) 


-29.1 


38.3 


±24.8 



TABLE X: The response of Si-0 bondlength and bond angles of a-cristobalite to r/n, 7733 and 
7723 strains. The number of different Si-0 bondlengths, Si-O-Si and 0-Si-O angles increase upon 
application of the rjn strain and further double under application of thhe r/23 strain. Angles are in 
degree and lengths in A. 

d/dr]ii d/dr]33 
(Si-Si) 0.190 -0.675 
Pj;(Si-Si) -0.226 -0.675 
P^(Si-Si) 1.676 

TABLE XI: The derivatives of the projection of the Si-Si vectors along the direction i (Pj (Si-Si)) 
with respect to strains rju and 7733 in a-cristobalite. 

with periodic boundary conditions used here. Its properties can then be studied theoreti- 
cally, also to get insight into the properties of amorphous silica. In fact, /9-cristobalite is the 
crystaUine phase of silica with density (2.18 g/cm^) and refractive index closest to those of 
amorphous silica. We have modeled /3-cristobalite in a cubic supercell containing 8 formula 
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units with M2d symmetry. We optimized the internal coordinates at the experimental den- 
sity of 2.18 g/cm^ (a = 7.13 A). A 4x4x4 MP mesh has been used in the BZ integration. 
The internal structural parameters of /5-cristobalite are assigned by the SiOSi angle, two 
OSiO angles and the Si-0 bond length. The calculated values are 142.9° (SiOSi), 107.5° and 
113.6° (OSiO) and 1.606 A (Si-0) in good agreement with the experimental values 146.7°, 
107.8° and 112.8°, and 1.611 A respectively^^. 

The dielectric tensor has two independent components eu = 622 and £33. The calculated 
dielectric constants and photoelastic coefficients are reported in table IXIII 

In the perspective to study the dielectric response of large models of amorphous silica, 
we have checked the convergence of the photoelastic coefficients as a function of the size of 
the mesh in the BZ integration for /5-cristobalite. Table IXTTl reports the results obtained by 
using the F-point only of the 24-atoms supercell and a mesh of 8 k-points in the first BZ 
which corresponds to the F-point only of a cubic supercell containing 192 atoms. The results 
with a mesh of 32 k-points are the "exact" values at convergence in the BZ integration. The 
dielectric and photoelastic tensor are already converged by using the F-point of the 192- 
atoms supercell. We have also computed the dielectric tensor and photoelastic coefficients 
of the 192-atoms supercell of /5-cristobalite within the Berry phase approach to electronic 
polarization, as developed by Putrino et ai.— . This latter calculation converges to the 
correct DFPT result, although with a different pace with the cell size by using the F-point 
only in the BZ integration, or, equivalently, with the fc-point mesb^. In fact, calculations 
within DFPT and BP with the same fc-points mesh give different results and the discrepancy 
decreases by increasing the fc-points mesh (or, equivalently, the cell dimension). In table 
IXIII the dielectric and photoelastic tensor of /?-cristobalite calculated within DFPT and BP 
are compared. The BP approach has been applied to a /3-cristobalite supercell containing 
192 atoms, by restricting the BZ integration to the F-point only, which is equivalent to the 
8 /c-points mesh of the 24 atoms supercell. However, the DFPT results with 8 k-points 
(24-atoms cell) are different from the BP results on an equivalent k-point mesh (F-point 
and 192-atoms cell). The BP approach convereges converges more slowly with cell size (or 
equivalently with the A;-points mesh) than DFPT. Moreover, the convereged value of e is 
approached from below within the BP approach. This effect is responsible for the difference 
in dielectric constant of the a-Si02 reported here (see infra) for a 81 supercell and a 162- 
atoms cell and for a 72-atoms cell recently considered in ref.— . 
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r 8 k 


32 k points 


BP 192 atoms 




2.337 2.281 


2.282 


2.181 


£33 


2.310 2.251 


2.252 


2.161 


P2I 


0.297 0.291 


0.292 


0.266 


P22 


0.131 0.124 


0.119 


0.165 


P23 


0.242 0.269 


0.267 


0.238 



TABLE XII: Dielectric and photoelastic tensors of /3-cristobalite calculated at convergence with 
respect to BZ integration (MP 4x4x4 mesh corresponding to 32 /c-points in the irreducible wedge) 
and with two other coarser meshes (P-point only and 8 A;-points in the IBZ, see text). The cubic 
unit cell containing 24 atoms has been used. The results obtained with the code CPMDi^ and the 
Berry-phase approach of reffi^ is reported for the 192-atoms supercell at the P-point only. 

From table IXIII it turns out that for system size of the order of 192 atoms the errors 
in the photoelastic coefficients due to the use of BP approach at P-point are of the same 
order of the errors introduced by the LDA itself. Since the BP approach is computationally 
less demanding that DFPT as implemented in the codes CPMD and PWSCF, respectively, 
we have used the former approach for the largest amorphous model described in the next 
section. 

V. AMORPHOUS SILICA 

We have generated several models of amorphous silica by quenching from the melt in 
classical MD simulations using the empirical potential by van Beest et alM as described in 
section II. We have considered simulation cells of two sizes: 81 and 162 atoms. We first 
report on the details of the quenching protocols and on the structural properties of the 
resulting amorphous models before discussing the dielectric properties of a-Si02. 

A. Structural properties 

The generation of amorphous models by quenching from the melt within Car-Parrinello 
simulations poses some unavoidable restrictions on the cell size and the quenching time, 
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which, in turn, negatively affect the quahty of the models produced. The main problem is 
that a very high quenching rates (in ref.— quenching rate was 10^^ K/s), prevents an adequate 
relaxation which freezes the concentration of small rings (three and four membered) to a 
value much higher than those expected in the experimental samples^. A further limitation 
is given by the small size of the systems that can be treated by ab-initio simulations (in 
ref.— the simulation supercell contains 72 atoms). 

In order to partially overcome these limitations, we have adopted a combined classical 
and ab-initio simulation scheme, that has been already reported to produce models of glassy 
Si02 with satisfactory structural and electronic properties, when compared to fully ab-initio 
models and to available experimental data^^. Within this approach, models of a-Si02 are 
generated by quenching from the melt by classical MD with the BKS potential at low 
quenching rate. The amorphous models thus obtained are then annealed by Car-Parrinello 
MD. Starting configurations of liquid silica have been obtained from unstable simple cubic 
crystal which rapidly reaches 8500 K. For both the 81- and 162-atoms models we have 
used an initially cubic box with density of 2.21 g/cm^. After 25 ps at 8500 K, the high 
temperature liquid is cooled to about 4000 K in 25 ps. The temperature at which atomic 
diffusion freezes over our simulation scale is much larger (~3500 K) than the experimental 
vitreous temperature (~1700 K). In order to obtain a good amorphous model, we have 
equilibrated the liquid at a temperature slightly above 3500 K and then quenched it slowly. 
The quenching protocol used to generate the 162-atoms model and a first 81-atoms model 
is shown in Fig. ^ It turns out that the quenching rate of Fig. ^ is still too fast to 
correctly reproduce the structural properties of amorphous silica in the small system (81 
atoms) which in fact appears to be strongly anisotropic. This problem is less severe for the 
larger system (162 atoms) which can escape from high energy local minima visited along the 
quenching. We have then generated another smaller model (81-atoms) with a much lower 
quenching rate by equilibrating the system for 5 ns at 3800 K and quenching to 300 K in 
2.5 ns (corresponding to a quenching rate of 1.4-10^^ K/s). The 81-atoms system referred 
to hereafter corresponds to this second quenching protocol. 

Ab-initio annealing for 1.0 ps at 600 K has been then performed by CPMD. For the 
smaller model, we have computed structural properties over a dynamical ab-initio micro- 
canonical run at 300 K, 0.6 ps long. The pair correlation functions and the bond angle 
distributions of the ab-initio and classical models at 300 K are compared in fig. El 
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FIG. 1: Quenching protocol adopted for the 162-atoms amorphous model. 




1 2 3 



r[A] 




angle [degree] 

FIG. 2: (a) Partial pair correlation functions and (b) bond angle distributions for the 81-atoms 
model of a-Si02, computed within classical MD (solid lines) and CPMD (dotted lines) at 300 K. 
Histograms in panel (b) refer to the structure optimized ab-initio at the theoretical equilibrium 
density and at zero temperature. 
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FIG. 3: (a) Partial pair correlation functions and (b) bond angle distributions for the 162-atoms 
model of a-Si02, computed within classical MD at 300 K. Histograms in panel (b) refer to the 
structure optimized ab-initio at the theoretical equilibrium density and at zero temperature. 
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23456789 10 23456789 10 11 
Silicon atoms per ring Silicon atoms per ring 

FIG. 4: The ring-size distributions of the 81-atoms (a) and of the 162-atoms silica models (b) 
computed according to Ref4^. 
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The average OSiO angle is very close to the tetrahedral value (109.5°), while the SiOSi 
bond angle distributions are quite different in the classical and CPMD simulations. The 
main effect of ab-initio annealing is a shift to lower angles of the SiOSi angle distribution. 
The Si-0 and 0-0 pair correlation functions are not affected, while the Si-Si pair correlation 
function reflects the modification in the SiOSi angle distribution, showing a longer tail at 
small Si-Si distances. The ab-initio annealing and optimization does not affect the topology 
acquired by the system during the classical MD quenching. The ring-size distribution of 
the two a-Si02 models is reported in Fig. IV Al The 81-atoms model produced with the low 
quenching rate has a low concentration of small rings (3- and 4-membered) with respect to 
the 162-atoms model and to previous models (72-atoms large) generated fully ab-initio^. 
After the ab-initio annealing performed with the softer Vanderbilt pseudopotentials, the 
final structure has been further optimized with norm- conserving pseudopotentials (for the 
calculations of the dielectric properties). The cell geometry has been optimized at fixed 
volume allowing orthorhombic distortions of the initially cubic supercell such as to produce 
a diagonal stress tensor. 

The residual anisotropy in the stress (a) is: 

^-422.8 -11.1 7.1 ^ 
(T= -11.1 -421.1 -2.6 (2) 
^ 7.1 -2.6 -422.2 y 

for the optimized ratios 6/a=1.051, c/a=1.002 in the 81-atoms supercell and 

/ _ /ion n 1 /I n q a \ 



a 



V 



-429.0 -14.0 8.4 
-14.0 -424.5 6.1 
8.4 6.1 -424.3 



(3) 



for the optimized 6/a=1.037, c/a=1.018 in the 162-atoms supercell. The large negative 
stress in eq. El is due to the so-called Pulay stress. The b/a and c/a ratio obtained in this 
way at the initial density of 2.21 g/cm^ have been then hold fixed and the volume varied 
to generate the equation of state reported in Fig. El As far as we know this is the first ab- 
initio calculation of the equation of state of a-Si02. The calculated E{V) points have been 
corrected for the discontinuities due to the incomplete basis set following the prescription 
given in ref.— and then fitted by a Murnaghan functional. The resulting equilibrium density 
{Peq), bulk modulus {B) and derivate of the bulk modulus with respect to pressure {B') for 
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-328.28 




Volume/atom [A^] 

FIG. 5: Ab-initio equation of state of the a-Si02 models, a) 81-atoms supercell. b) 162-atoms 
super cell. 

the 162-atoms and 81-atoms models are peq = 2.34 g/cm^, B = 39.0 GPa, B' = —7.01 and 
Peg = 2.31 g/cm^, B = 46.2 GPa, B' = —6.28, respectively. These results, especially for 
the larger model, are in good agreement with the experimental values of peg = 2.205 g/cm^, 
B = 36.7 GPa, B' = —5.31—. Note in particular that the highly unusual negative value of 
B' is well reproduced by our calculations. The bulk modulus would probably improve by 
optimizing b/a and c/a at each volume, held fixed in the equation of state of Fig. The 
error in B is smaller for the larger system which has less geometrical constraints and can 
more easily respond to the density change. In fact, we have also computed the equation of 



size (atoms) 


81 162 648 


Peg (g/cm^) 


2.23 2.30 2.29 


B (GPa) 


53.5 40.5 34.1 


B' 


-9.0 -4.3 -6.3 



TABLE XIII: Equilibrium density (peg), bulk modulus (B) and first order derivative of the bulk 
modulus (B') of amorphous silica models optimized with the BKS potential. 

state within the classical molecular dynamics scheme, using the empirical potential of ref.— 
for the 81-, 162-atoms models and for a larger amorphous model containing 648 atoms. The 
results, reported in table IXTTH indicate that the overestimation of the bulk modulus in the 
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ab-initio calculations can be ascribed mainly to a size effect. 

Bond angle distributions of the fully optimized system at the equilibrium density are 
plotted in figures |21 and El for the 81-atoms and for the 162-atoms model, respectively. The 
shape of the SiOSi angles distribution strongly depends on the quenching rate. The model 
quenched at lower rate (81-atoms cell) has a larger mean SiOSi angle both in the classical 
simulation and after ab-initio annealing. Furthermore, for this model the concentration of 
small (<120'') SiOSi angles is lower than in the model generated fully ab-initio in ref.— and 
fits better to experimental data. 



B. Dielectric and photoelastic properties 



Since the photoelastic coefficients are very sensitive to density, we have computed the 
dielectric properties of amorphous silica at the theoretical equilibrium volume (cfr. previous 
section). The residual, small anisotropy in the structure can also be identified by inspection 
of the dielectric (e) tensor: 

^ 2.297 0.013 -0.009 ^ 

0.013 2.292 0.007 , (4) 
\ -0.009 0.007 2.288 j 
for the 81-atoms supercell and 

^ 2.192 0.009 -0.004^ 

0.009 2.198 -0.018 (5) 
\ -0.004 -0.018 2.191 J 
for the 162-atoms supercell. For a fully isotropic and homogeneous system, as expected for 
a truly amorphous system, e should obviously reduce to a single number. The experimental 
value of e at p =2.2 g/cm^ is 2.125^. At the same density the theoretical dielectric constant 
(1/3 Tr(e)) is £=2.260 for the 81-atoms cell and e=2.148 for the 162-atoms cell. The misfit 
with respect to experiments is similar to those already found for the crystalline phases of 
silica. For the largest system (162-atoms) the use of F point only within the BP framework 
partially compensates for the error due to LDA, as occurs for /3-cristobalite (cfr. section 
3.5.3) resulting in an accidental better agreement with experiments. The dielectric constant 
computed at F within the BP for a 72-atoms cell of a-Si02 in ref.— is 2.00, which confirms 
that, within BP, e converges from below by increasing the cell size. 
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81-atoms 






162-atoms 




d/drjii 


d/dri22 




d/ dr]22 


d/dr]3s 


Si-0 


0.151 


0.154 


0.157 


0.15 


0.11 


Si - 6 - Si 


68.6 


69.1 


73.4 


68.5 


57.5 


Px (Si-Si) 


1.488 


-0.037 


-0.122 


-0.105 


-0.053 


P,(Si-Si) 


-0.141 


1.553 


-0.049 


1.452 


-0.044 


P. (Si-Si) 


-0.124 


-0.067 


1.490 


-0.098 


1.355 



TABLE XIV: The derivatives of the SiO distance, SiOSi angles and Si-Si vector distances with 
respect to strains r/n, 7722 and 7733 in the models of a-Si02. Pj (Si-Si) denotes the projection of the 
Si-Si vector distance on the i-th axis. 

At the theoretical equilibrium density the dielectric constant is e=2.292 and e=2.194 for 
the 81- and 162-atoms cells, respectively 

We have computed the response to strain of the structural properties of the amorphous 
models as reported in table IXIVl 

The calculated photoelastic coefficients of the 81- and 162-atoms supercells at the theo- 
retical equilibrium density are reported in table IX Vl The photoelastic coefficients have been 
calculated by finite differences from the dielectric tensor of systems with strains from -1 % 
to +1 %. As discussed in section II and IVb the calculation on the 81-atoms cell have been 
performed within standard DFPT, while the calculations on the larger 162-atoms supercell 
have been performed within the Berry phase. Residual structural anisotropics are responsi- 
ble for the differences among pu, P22 and P33 values and among pu, P21 and P32. A measure 
of the model isotropy is given by the fulfillment of the Cauchy relation ( 2^44=^11 —pu)- 
The value of P44 computed by applying the 7723 shear strain in the 81-atoms model is -0.074 
(cfr. tab. IXV|) which compares well with the value of {pu —pi2)/2, ranging from -0.069 
to -0.092 which suggests a sufficient homogeneity of the small model as well. Agreement 
with experiments is fair, but less satisfactory than the performances of DFPT-LDA for crys- 
talline phases. In particular, the ratio pn/pu is largely underestimated in our model with 
respect to experiments. In view of the excellent results of DFPT-LDA for crystalline silica, 
the discrepancy with experiments for a-Si02 could be naturally ascribed to an insufficient 
quahty of our a-Si02 models, either due to a size effect or to a still too fast quenching rate. 
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81-atoms 162-atoms exp. 





0.072 




0.125" - 0.100^ 




0.045 


0.047 




Fid 


0.053 


0.085 




7)1 o 


0.217 


0.223 


27" - 285^ 


T)oi 


0.230 








0.224 






P32 


0.209 


0.227 




Pl3 


0.209 


0.238 




P23 


0.197 


0.235 




P44 


-0.074 




-0.073" 



TABLE XV: Photoelastic coefficients of the amorphous sihca models containing 81 and 162 atoms 
compared to experimental data, a: ref.-^, b: ref4i. 

However, we should also note that the measurement of the absolute values of pn and pu 
is very delicate and to our knowledgment only two independent measurements have been 
performed so fai*^^^. 

In order to obtain insight into the different microscopic contributions to photoelasticity 
and shed light on the relationship between structural properties and photoelastic response, 
we have developed a phenomenological model of photoelasticity in silica, as described in the 
next section. 

VI. PHENOMENOLOGICAL MODEL OF PHOTOELASTICITY IN SILICA 

In order to identify the main microscopic features ruling photoelasticity in silica, we have 
devised a phenomenological model of the dielectric properties transferable both to the crys- 
talline and amorphous phases of Si02 by fitting on the ab-initio database on photoelasticity 
in crystalline silica reported in the previous sections. 

At this aim, we have assumed that the dielectric response of silica could be embodied 
in an ionic polarizability tensor of the oxygen ions, whose value is assumed to depend on 
the SiOSi angle only. The dielectric tensor e can be obtained from a site dependent oxygen 
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polarizability, Oi, as described below. The dipole moment pi at site i is related to the local 
field Ei via the polarizability tensor a as 

Pi = a^Ei . (6) 

The local electric field E^ acting on a polarizable object at site i is then^: 

N 

E, = D- Y.Tl-pk, (7) 

k^i,R 

where D is the electric displacement, the total number of polarizable centers per cell and 
T-f is a symmetric tensor defined as: 



Ti = ViVfc 1 — I = -3- 



1-3^ 



(8) 



where R are the Bravais lattice vectors and rik is the distance between sites i and k in cells 
separated by R. The sum of T^f over the Bravais lattice vectors R is performed as an Ewald 
sum as described in ref.— . 

The electric displacement D is in turn equal to D = E + AttP, where E is the macroscopic 
electric field E and the electric polarization P can be expressed as 

1 AT 1 AT 

P=yT.P^=yT.2l^^^ (9) 

i i 

where V is the cell volume. A set of linear equations which relates the local field Ei to the 
macroscopic electric field is obtained from eqs. [3121 as: 

E = E,-j:(y-j:Ti\a^E„ (10) 

^y^^ \ R / 

or in a compact matrix form 

E = a-Rm} . (11) 

/ is the identity matrix, E and Ei are expressed as 3N vectors and 5 is a 3Nx3N matrix 
consisting of 3x3 blocks Bjk defined as 

For a cubic lattice equations El yield the Clausius-Mossotti (Lorentz-Lorenz) formula (see 
ref.-^ for details). However, none of the crystalline polymorphs of Si02 is cubic, thus the 
term containing the dipole lattice sum (X) T^^) has to be explicitly evaluated. 
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From eq. IHl and eq. El the following relation is obtained: 

p = I^T.^{l-R).'e. (13) 

The inversion of the matrix — is performed in the 3N-dimensional space. Finally, the 
dielectric susceptibility tensor x (P = xE) results in: 

i-yj:^{L-R):. (14) 

From X, both the dielectric {e = 1 + 47rx) and the photoelastic tensors can be derived. 

In order to exploit eq. El in our empirical model we have assumed that the polarizable 
centers in silica are the oxygen ions. Secondly, we have also assumed that the oxygen polar- 
izability depends on the SiOSi angle only. These choices can be justified only a posteriori 
from comparison of the results produced by the phenomenological model with the ab-initio 
data. However, to support our model we can argue that most of the valence charge in Si02 
is located on the oxygen ions and, secondly, that both the structural differences among the 
different tetrahedral phases of silica (either crystalline or amorphous) and the deformations 
induced by strain mainly rely on the distribution of the SiOSi angles. The internal structure 
of the tetrahedra (OSiO angles and Si-0 bond lengths) is almost the same in the different 
polymorphs and undergoes minor deformations upon strain since the open structure of silica 
allows to accommodate strain mostly by rotations of the corner-sharing tetrahedra. 

Our model of silica polarizability can also be understood as a particular case of the 
commonly used bond polarizability model (BPM). In the BPM the polarizability is expressed 
as the sum of Si-0 bond polarizabilities given by^ 

1 / T 'T ' \ \ 

o^ij = 3^"^ + 2aT)^ii + ("i ~ "r) ( ^ - -Sij] (15) 

where f = fsi—ro is the vector joining the silicon and oxygen atoms, and and represent 
a longitudinal and a transversal polarizability of the Si-0 bond, ai and ar usually depend 
on the Si-0 bondlength (see fig. IH)). In our model we have introduced a second transverse 
bond polarizability, ot', in the direction orthogonal to the Si-O-Si plane (direction z in fig. 
inj. We have made a/,, ax and a^' dependent on the SiOSi angle {6), but independent 
on the Si-0 bondlength. The dependence of the polarizability on the Si-0 distance is 
crucial to reproduce Raman intensities, which are ruled by the modulation of the bond 
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FIG. 6: Sketch of the Si-O-Si unit and of the bond contributions to polarizabiUty. 



a 



polarizability upon phonon displacements, that obviously involve both angle modulations 
and bond stretching. Conversely, the change of a/,, ay and aT> with Si-0 bondlength is 
not essential to describe photoelasticity, since strain-induced structural relaxations mainly 
consist of bond angles deformations not affecting the Si-0 distances. 

In this approximation the dielectric response can be cast into a model based on the 
polarizability of the oxygen ions only which is given by 

/c(^) +7cos2(^/2) ^ 

c(^) +7sm2(^/2) (16) 

y oit\Q) I 

with c = 1{a.L + ot + otO/S and 7 = a/, — «t (cfr. fig. IHl for axis orientations). The 
contribution of each polarizable unit to the dielectric susceptibility is 

= ^a{0)Ri (17) 

where B4 is the rotation matrix that operates the transformation from the local reference 
system represented in figure El to the absolute reference system of the solid, in which the 
i-th Si-O-Si unit is embedded. The parameters c, 7 and a^' have to be determined and, in 
principle, all of them depend on 6. We have decided to fit these parameters on the dielectric 
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tensor of a-cristobalite at several densities (see fig. [7j), since a-cristobalite responds to com- 
pression mainly modifying the Si-O-Si angles, maintaining the Si-0 bondlength unchanged^. 
We have checked that changing the density of a-cristobalite from 2.83 to 2.05 g/cw? (the 
theoretical equilibrium density is 2.54 g/cm^) the Si-O-Si angle changes from 130° to 165° 
with a maximum variation of the Si-0 bondlength of 0.006 A (see fig. |H1)- 
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FIG. 7: The dielectric tensor of a-cristobalite as a function of the Si-O-Si angle. 
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lattice parameter [a.u.] 

FIG. 8: The modulation of the Si-O-Si angle in a-cristobalite as a function of the lattice parameter 
a. 

However, for a-cristobalite the equations for x (sq- IT^ independent by symmetry are 
only two, not sufficient to fit all the three functions c, 7 and ut' in eq. ^1 For this reason, 
7 has been assumed independent on 9 and set equal to the value obtained from the fitting 
of the ab-initio Raman spectrum of a-quartz within the BPM by Umari et ai.— (7 = 9.86 
a.u.). The parameterization of the BPM proposed by Umari et ai. reproduces the ab-initio 
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Raman intensities of a-quartz within 15%. The results for c{6) and ax'iO) obtained from 
eq. 1131 have been interpolated by a second order polynomial (see fig. inj, whose coefficients 
are reported in table IVll 





ao ai 


a2 


c{e) 


-20.46 18.38 


-2.942 




-6.141 9.956 


-1.730 



TABLE XVI: Coefficients of the polynomials that fit c{9) and aT'{d) {0 is expressed in radiants). 




FIG. 9: The functions c{9) (solid line and circles) and aT'{0) (dashed line and squares), which 
assign the oxygen polarizability (see text). The data are the result of the fitting on the dielectric 
tensor of a-cristobalite at diff'erent densities. The axx and ayy components of the polarizability 
tensor (dashed lines) as obtained from eq. 1161 

As shown in fig. IHlthe oxygen polarizability increases by increasing the Si-O-Si angle. In 
fact, the Si-0 bond becomes more ionic by increasing the Si-O-Si angle and the larger charge 
on oxygen ion makes it more polarizable. The increased ionicity of Si-0 bond at larger SiOSi 
angles is further confirmed by the dependence of the effective charge of oxygen ions on the 
SiOSi angles shown in^J The latter figure reports the distribution of the effective charges 
(l/3TrZ*) of oxygen ions in different sites of our a-Si02 model, computed within DFPT— . 

The transferability of the oxygen polarizability shown in fig. |H1 has been checked by 
comparing the dielectric constants of the other silica polymorphs discussed in the previous 
sections, as obtained within the phenomenologial model and by DFPT (see table lyTj) . The 
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FIG. 10: Dependence of the Born effective cliarges of oxygen ions on the Si-O-Si angle for our 
81-atoms model of a-Si02. 
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2.694 
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£11 
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P=7 GPa 
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2.865 
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£11 


2.281 


2.252 
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2.251 


2.227 




£11 


2.297 


2.285 


a-SiOa 


£22 


2.292 


2.281 




£33 


2.288 


2.275 



TABLE XVII: Comparison of the dielectric constants computed within ab-initio DFPT and within 
the phenomenological model fitted on a-cristobalite. 

agreement is overall satisfactory. In particular, birefringence of quartz and /?-cristobalite 
are well reproduced. The discrepancies for quartz at high pressure may be due to the 
change of the Si-0 bondlength upon compression (see table |V|), whose effect on the oxygen 
polarizability has been neglected in our model. This effect is less important in a-cristobalite, 
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whose more open structure can easily accommodate density changes with minor effect on 
the Si-0 bondlength. The phenomenological model reproduces the dielectric constant of 
a-Si02 as well. 

The transferability of the model has been further checked by computing the photoelastic 
tensor. This is given by the derivative of the dielectric susceptibility tensor x (gQ- fT^ with 
respect to the strain tensor rixpL as 



brackets indicate 3x3N matrices. The change of the oxygen polarizability with strain can 
be expressed in terms of dc{9)/d6 and daT'{0) / 06 deduced from Fig. 121 AH the other 
derivatives have been obtained by finite differences. Eq. ^Jcan provide an estimate of the 
photoelastic coefficients of all silica polymorphs made of corner-sharing tetrahedra. 

The photoelastic coefficients calculated within the phenomenological model for a few 
silica polymorphs are compared to ab-initio data in table IXVlllI (second and first column, 
respectively). The same comments made above for the dielectric constants still hold for 
the photoelastic tensor. The agreement with ab-initio data is better for system at ambient 
conditions than at high pressure. Overall we can conclude that the phenomenological model 
has good transferability. 

Table IXVlllI reports also the contribution from the individual terms supplying to the 
expression of dx/drjx^ in eq. [THl 

The first term in the right-hand side of eq. ^1 (— X(5a^) comes from the derivative of the 
density, which is null for strains preserving the volume i.e. for X ^ fi. The contribution 
of this term is always positive. The second term in eq. ^1 contains the derivative of the 
polarizability of the Si-O-Si units. It can be expressed as a function of the Si-O-Si bond 
angle and depends on the derivative of the functions c{9) and aT'{0): 




(18) 



where the matrix A is defined as A = (/ — 5) ^ (see eq. IT^ . and the arguments in square 



da{9) _ daje) dO 
dr] 09 Or] 



(19) 



32 



DFPT model dR/dr] Lorentz da/ 89 = 
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P=0 GPa 
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a-quartz 
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0.239 
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0.049 


0.290 


/3-cristobalite p2i 
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0.275 


0.040 


0.236 


0.470 




P31 


0.269 


0.236 


0.041 


0.207 


0.351 




Pll 


0.072 


0.103 


-0.082 


0.074 


0.168 


a-Si02 


P21 


0.230 


0.217 


0.042 


0.227 


0.319 




P31 


0.224 


0.214 


0.038 


0.221 


0.306 




Pll 


0.218 


0.161 


-0.013 


0.100 


0.353 




P21 


0.244 


0.276 


-0.010 


0.275 


0.363 


a-cristobalit 


e P31 


0.293 


0.319 


0.025 


0.289 


0.408 




Pl3 


0.293 


0.325 


0.046 


0.296 


0.449 




P33 


0.152 


0.154 


-0.091 


0.105 


0.228 



TABLE XVIII: Comparison of the photoelastic coefficients yielded by the phenomenological model 
and by ab-initio DFPT. The data on a-Si02 refers to the 81-atoms model. See text for the 
description of the other columns. 



da{9) 
86 



\ 






daj,,{e) 



(20) 



with 

( 0.5^stn{e) n \ 



This term provides a negative contribution to all the components of the photoelastic tensor. 
The last column in table IXVIIII reports the photoelastic coefficients obtained from eq. 
but neglecting the contribution from da/ 89. By comparing the second and the last columns 
in table IXVllll it is clear that the change in polarizability with the SiOSi angle is essential 
to correctly reproduce the ab-initio data. Although it has not been explicitly displayed in 
eq. UHl the derivative of a also contributes to the term dA/drjx^. 
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The last term in eq. ITHl depends on the derivative of the rotation matrices Ri with respect 
to strain. It is related to the geometrical effect of alignment of the Si-Si axis of the Si-O-Si 
unit (cfr. fig. ^ along the axis of a tensile strain. This effect results in an increase of the off- 
diagonal components and a reduction of the diagonal components of the photoelastic tensor. 
This effect accounts for most of the observed difference between diagonal and off-diagonal 
photoelastic coefficients (pn and pi2 for instance). This is proven by the results in the third 
column of table IXVIIIl which reports the contribution to the photoelastic coefficients due to 
the last term in eq. ITHlonlv. In fact the difference pu —pu for a-Si02 and /?-cristobalite, and 
P33—P13 for a-cristobalite and a-quartz, in the second column of table IX VI I H are close to the 
results in the third column. We note that the identification of the alignment of the Si-O-Si 
units as the source of a large value of pn —pu and P33 —pis is consistent with the structural 
response data in tables IVIIl IXII and IXIVI The small value of pu — pi2 in a-cristobalite is 
consistent with the weak e alignment of the tetrahedra since for strain along the x-axis is 
weak (cfr. table IXT|l . 

The term depending on the derivative of the matrix A accounts for the the presence of 
local fields correction and its modification upon strain. Its contribution to the photoelastic 
coefficients is inferred by computing pij within the Lorenz-Lorentz approximation. The 
results, shown in the 4th column of tab. IX VI I H demonstrate that local field effects beyond 
the Lorenz-Lorent approximation (i.e. inclusion of anisotropic local fields produced by near 
molecules) modify the photoelastic coefficients for crystalline phases up to 45%. As expected, 
corrections to the Lorenz-Lorentz approximation are small for the amorphous silica models. 
The phenomenological model allowed us to address the role of the size of the simulation cell 
on the photoelastic tensor of a-Si02. This has been done by computing the photoelastic 
coefficient within the phenomenological model for the same 81-atoms and 162-atoms models 
discussed so far and for a larger model of 648 atoms, all optimized with the BKS potential. 
The results, reported in table IXIXl quantify the errors due to finite size effects. Isotropy in 
the diagonal components pu are achieved only for the largest model. However, errors due 
to finite size effects for cells containing 81-162 atoms are only slightly larger than the errors 
due to LDA-DFT. 
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81-atoms 162-atoms 648-atoms 



i^ll 


0.115 


0.100 


0.132 


P22 


0.161 


0.143 


0.141 


P2I 


0.227 


0.231 


0.234 


P31 


0.237 


0.237 


0.227 


P12 


0.245 


0.216 


0.220 


P32 


0.240 


0.239 


0.220 



TABLE XIX: Photoelastic coefficients computed within the phenomenological model for three 
models of a-Si02 optimized with the BKS classical potential. 

VII. DISCUSSION AND CONCLUSIONS 

Photoelasticity of crystalline and amorphous silica has been studied by density func- 
tional perturbation theory. This framework has been checked successfully first on paradig- 
matic semiconducting and insulating systems, silicon and MgO. The same framework has 
been applied to a-quartz providing photoelastic constants in good agreement with previous 
calculations^ and with experiments. The photoelastic coefficients have been computed also 
for a-quartz under pressure up to 7 GPa, and for a-cristobalite and /5-cristobalite at normal 
conditions. The theoretical results obtained for these latter systems are predictions open to 
experimental verification. 

We have thus verified that the simple DFPT-LDA framework is suitable to reproduce 
accurately the photoelastic coefficients. Correction of the DFT-LDA band gap by a scissor 
operator is necessary to reproduce the dielectric constants but has minor effects on the 
photoelastic tensor. 

We have then studied photoelasticity for models of amorphous silica, containing up to 
162 atoms, generated by quenching from the melt in molecular dynamics simulations. The 
calculated photoelastic coefficients are in good agreement with experimental data. Our 
ab-initio framework has thus valuable predictive power also for the amorphous phase. 

Along with the calculated photoelastic coefficients, we have collected data on the struc- 
tural response to strain of different crystalline phases. The analysis of these data has aided 
us in the development of a phenomenological model of photoelasticity in silica polymorphs. 
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In fact the analysis of the structural changes upon strain has shown that silica responds to 
strain mainly via rotations of the tetrahedral units (thus changing the Si-O-Si angles) with 
smaller changes of the intra-tetrahedral (0-Si-O) angles. For instance, the most prominent 
effect of tensile stress is the alignment of the Si-O-Si units to the direction parallel to the 
strain axis. 

These results suggested us to model the dielectric properties of silica polymorphs in 
terms of Si-O-Si polarizable units with polarizability dependent on the Si-O-Si angle only. 
At this aim we have developed a phenomenological model fitted on the dielectric constant 
of a-cristobalite at different densities. The model consists of polarizable oxygen ions with 
anisotropic polarizability dependent on the Si-6-Si angle. The validity of the model has 
been proven by comparing dielectric and photoelastic coefficients computed within the phe- 
nomenological model to the results of DFPT calculations. The agreement is good and the 
simple phenomenological model sheds light on the main contributions to photoelasticity in 
silica polymorphs. It turns out that the change in the oxygen polarizability with strain, via 
the SiOSi angle is essential to reproduce the photoelastic response. The neglect of this term 
introduces errors of up to 50 %. Secondly, the difference in the diagonal and off-diagonal 
photoelastic coefficients (e.g. pu and P21), which is particularly large in a-Si02, is due to 
the alignment of the SiOSi units along the axis of tensile strain. Finally, while local fields 
beyond the Lorentz-Lorent approximation (via dipole-dipole interactions) and their change 
upon to strain are important to reproduce the dielectric and photoelastic properties (expe- 
cially birefringence) of the low-symmetry crystalline phases, they are obviously negligible for 
a-Si02. This result confirms that our models of a-Si02, despite their small size, are already 
good model of a homogeneous isotropic amorphous network. 

In summary, our work has demonstrated the reliability of the ab-initio DFT-LDA frame- 
work in the study of photoelasticity of silica systems and has provided a phenomenological 
model of the dielectric properties transferable to several crystalline and amorphous poly- 
morphs. The same ab-initio framework, supplemented by the phenomelogial model, could 
also provide a microscopic description of the photoelasic properties of more complex silica- 
based glasses. 
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